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Abstract. We propose a new application of the scheme of fiber bundles to the 
counting statistics. Using the fiber bundle discussions, a cumulant generating function 
for current in stochastic processes is divided into two parts, i.e., the dynamical 
quantity and the geometric one. We will show that the introduced geometric quantity 
is connected to a kind of excess quantities, which have been investigated a lot in 
nonequilibrium physics. Numerical experiments for a specific nonequilibrium model 
are performed, and we confirm that the geometric quantity has the characteristics of 
the excess quantities. 
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1. Introduction 

Studies of nonequilibrium systems liave been focused a lot mainly in the last decade. One 
of the main topics is the transitions between stationary states. During the transitions, 
a system causes a kind of flows or currents of physical quantities. For example, as for 
a transition between two equilibrium states, the change of thermodynamic quantities 
can be calculated using equilibrium thermodynamics. However, different from the 
equilibrium systems, the nonequilibrium systems have steady flows or currents, which 
stem from their nonequilibrium nature, even if we do not consider any transitions 
between stationary states. In order to investigate such nonequilibrium systems, a 
concept of 'excess' quantities plays an important role [THE]. That is, if we change 
the system using external forces, the nonequilibrium systems can cause extra currents 
caused by the change of states, in addition to the contributions from the steady currents. 
For such excess quantities, there are various studies, but the definition of the excess 
quantities is not obvious and there are many proposals for the excess quantities. 

Recently, the geometric phase concept has been used to investigate flows or currents 
in stochastic processes. For cyclic perturbations, the geometric phase concept has given 
successful results P-E]. However, in order to study the excess quantities for general 
cases, it is needed to treat noncyclic schemes; for example, we consider a sudden change 
of the transition rates of a stochastic system, and then the system does not return 
to the initial state. When the change of the transition rates is very slow (so-called 
adiabatic cases), a path-integral formulation is available and a quantity called a noncyclic 
geometric phase has been introduced. To our knowledge, at least two different definitions 
have been proposed [T5l[T6]. As for the nonadiabatic definition of the geometric 

phase has been introduced from a mathematical viewpoint [17]. However, it has been 
unclear whether the proposed geometric phase in [17] has physical meanings or not. 

In the present paper, a new geometric quantity is proposed, based on the scheme 
of fiber bundles. Although the scheme is similar to the previous work [17], we will show 
that the proposed quantity in the present paper has the characteristics of the excess 
quantity, and hence it would be meaningful physically. 

The present paper is constructed as follows. In section 2 we explain the problem 
settings of the counting statistics. Section 3 gives the application of the scheme of fiber 
bundles, and a geometric quantity is introduced. In section 4, we perform numerical 
experiments for a nonequilibrium model; the numerical results confirm that the proposed 
geometric quantity has the characteristics of the excess quantity. Section 5 gives some 
concluding remarks. 

2. Problem settings 

We consider a stochastic process with finite discrete states; the number of states is S. 
The dynamics of the stochastic process is described by the following master equation: 

^Mt)) = K{tMt)), (1) 
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where \p(t)) = {pi{t) , p2{t) , . . . ,ps(t))'^ is a vector for the probabihty distribution of the 
system and K{t) = {Kss'{t)} the transition matrix; the element Kssi{t) is related to the 
transition rate for the state change s' — s at time t. For generality, we consider that the 
state change s' — ?■ s is caused through several ways (see the example in Section 4). For 
each state change s' — )■ s, we consider the transition rates w\g,{t),w1^i{t), . . . yW^f {t); 
totally Vss' ways of transitions exist for the state change s' — s. In addition, we set 
Vss = 1 for all s and wl^it) = 0. Using the transition rates {'w'^g,{t)}, the elements of 
the transition matrix K(t) is written as follows: 

K,At)^Y.w:At)-\sm,s', (2) 
i/=i 

where 6s is the Kronecker's delta and 



s'=l u=l 

We next consider 'currents' in the stochastic process. One example of the currents 
is the number of the occurrence of specific transitions (for example, see [9]), and another 
one is the entropy production (for example, see [IE]). Here, we consider a quantity a^^, 
defined for each state change s' ^ s and the way of the transition; if we have one 
state change s' -> s by the i^-th transition way, we count a quantity a'^^,. Our aim 
is to calculate the statistics of an accumulated quantity. In order to formulate the 
accumulated quantity, we fix a time interval [0, T] and consider a trajectory of the state 
change, [s{t)]f^Q, where s{t) G S. The trajectory is characterized by jump events at 
discrete times (^1,^2, • • • ^^n-i), and the trajectory [s(t)]^Q is expressed as s(t) = Si for 
ti~i < t < ti with to = and t„ = T. In addition, the jump event at time ti is specified 
by the z/j-th way of transitions for the state change s' s. Then, the accumulated 
quantity is defined as 

n-l 

where z/j G {1, . . . , Ki+isi}- The generating function for the accumulated quantity A{T) 
is given by 

F(x,T) = (e'^^(^)>, (5) 

where the average (■) is taken over realizations of trajectories with the initial distribution 
\p(t)). Derivatives of F{xyT) with respect to x give various information about the 
statistics of the accumulated quantity. 

The problem here is: How should we calculate ([5])? One may consider that it is 
difficult to take the average over the trajectories. In order to evaluate ([5]), the scheme of 
the counting statistics has been developed well and used in various context (for example, 
see [IE])- The consequences of the counting statistics give us the following simple way 
to evaluate the generating function. For details of the derivations, see [19], for example. 
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We here consider the following matrix; 

HsAx,t) = ^i^;,,(t)exp(x<,0 - Xsm,s'. (6) 

Note that the matrix H{x,t) = {Hss'ix^t)} is not a 'probability' transition matrix. 
Using the matrix H{x,t), the following time-evolution equation is considered: 

^^mx,t)) = Hix,tmx,t)). (7) 

Note that the vector is not related to the probability distribution; using a 

vector (1| = (1, 1, ... , 1), we have 7^ 1 in general. 

The remarkable consequence of the counting statistics is as follows: Using the 
solution of ([7]) at time T with an initial condition \ip{Xi't = 0)) = \p{0)), the generating 
function in (|5]) is given by 

F{x,T) = {mx,T)). (8) 



3. Fiber bundles for counting statistics 

All statistics of A{T) defined by (jl]) are given by the generating function ([8]), and here 
we rewrite the scheme by using the fiber bundle discussion [20||2T|. Notice that although 
the basic scheme is similar to [IT], the following definitions and discussions adequately 
give a quantity with the excess property. 

Firstly, we introduce a 'bra' vector {iplx, t) \ and define the time-evolution as follows: 

^^{^ix,t)\ = mx.t)\H{x,t), (9) 

where the initial condition is {ipix^t = 0)| = Note that —H{x,t) has been used 
to define the time-evolution in [17], and the difference of the time-evolution operator is 
important to obtain a quantity which seems to be meaningful physically. Additionally, 
{■ilj{Xyt)\ ili^iXjt)))^ in general because of the non-Hermite property of H{x,t)- (Of 
course, the initial condition is also different from that of \ip{x, t)).) Using the 'bra' vector 
{■ip{Xyt)\ ^-nd the 'kef vector t)), the generating function F{x^T) is expressed as 

F(x,r) = (^(x,o)|^(x,T)). (10) 

Secondly, we introduce new state vectors which is useful to define the fiber bundles. 
Note that the inner product of the bra vector (■?/'(x,t)| and the ket vector |'?/'(x,t)) is 
not conserved in general because ^ ((^/^(x, t)|'?/'(x, t))) = 2(-?/'(x, t)|iJ|-?/'(x, t)). Instead 
of (^^(x, ^)| and |-?/'(x,^)), we introduce the corresponding vectors whose inner product 
is explicitly conserved. That is, defining a quantity 5(X;^) as 
.,._ {i^{x,t')\H{x,t'mx,t')) 

^""^ mx,t')mx,t')) ' ^''^ 

the following new bra and ket vectors are introduced: 

|0(x,t)) = exp(^-^*5(x,t')dt') mx,t)), (12) 
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=exp 6{x.t')dt'j (13) 

Hence, the time-evolutions of these vectors are given as 

^|0(x,t)) = [H{x,t) - Six,t)] mx,t)), (14) 

^^{^ix,t)\ = {vix,t)\mx,t)-6ix,t)]. (15) 
Note that the inner product is conserved because 

^(Mx,t)l0(x,t))) = o, (16) 

and from the initial conditions, we have 

(^(X,t)|0(x,t)) = l. (17) 

Thirdly, we consider the principal fiber bundle P whose fiber F is given by 
M+ = {x E M|x > 0}. The base space M is constructed as follows: P is constructed 
by the set of all states {|0(x, ^))}, and we introduce the equivalent relation ~ as a kind 
of proportional relation; if |0(x,t)) = c|0'(x,i')) with c > 0, |0(x,t)) and \(p'{Xyt')) are 
considered as equivalent ones. We therefore have the base space M = P/ ~. Using a 
natural projection map tt : P — t- M, the construction of the principal fiber bundle is 
finished. We also consider its dual constructed by {(v9(x, t)|}. 

Fourthly, the following connection is introduced for P. Let \4>{x, s)) be a curve in 
P, and write the tangent vector to this curve as \u) = ^|0(X)'5))- Then, using a dual 
curve (</)(x,s)|, the connection is defined as 

^(^^^^^M^^ifMM). (18) 

For the dual bundle, we introduce a slightly different connection as follows: 

(<^(x,s)|0(x,s))' 

where {v\ = -j^ifix^ Note that (v9(x, s)|0(x, s)) 7^ 1 in general. In the present 
paper, we treat only cases with conserving inner products; (</5(x, ■s)|0(x, s)) is constant. 
That is, when we have a curve |(/)(x, s)) in P, the dual curve (v9(x, s)| is selected as 
^(((^(x, s)|0(x, s))) = 0. Actually, the time-evolutions in (fT4|) and ( !T5|) conserve the 
inner product, and geodesic curves discussed later are assumed to have this property. 
Hence, we have A{x,s) = y4^(x, s). Additionally, in the time-evolution of (HM . the 
connection is always zero because 

Mx,t)l^l0(x,t)) = O. (20) 
Furthermore, the gauge transformation is given as 

|0(X,3)) ^e^«|0(x,s)), {^{x,s)\ ^ {^{x,s)\e-^^'\ 
A{x,s)^A{x,s) + ^, (21) 
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where /3{s) e R. 

Why do we have to use the above settings? The reason is that using the scheme of 
the fiber bundle, we can naturally define a 'geometric' quantity which are determined 
only by the shape of the trajectory of the time-evolution. That is, for a closed curve in 
P, the following geometric quantity is introduced: 

^{X,T) = j A{x,s)ds. (22) 

Although one may consider that the integral defined by the closed curve is meaningful 
only when we consider a cyclic evolution for (fT4|) . the following construction makes the 
integral (|22l) possible even for noncyclic time-evolution cases. The initial vector 0)) 
and the final vector T)) does not make a closed curve in general cases; there is only 
a one-way path from |</)(x, 0)) to \(f){XiT)) obeying the time-evolution f|T^ . In order to 
make a closed curve from the one-way path, a geodesic curve is useful. The basic scheme 
has been written in [I71[22] , and we explain the details of the usage of the geodesic curve 
in Appendix. The important point here is that the integral on the closed curve is given 
as follows: 

(fA{x,s)ds=[ A{x,s)ds+ [ A{x,s')ds' (23) 

J J (time-evolution) J (geodesic curve) 

= ln(^(x,O)|0(x,T)). (24) 

The first term in (|23l) vanishes because the connection defined by ( fTSl) is zero on 
the time-evolution in iHM and ( IT5l) : only the second term contributes and gives 
Although A{x, s) in the first term in f l23|) becomes non-zero when one considers a gauge 
transformation, the geometric quantity 'y{x,T) is gauge invariant. 

Using the introduced quantities S^Xyt) and 7(x, T), the generating function can be 
rewritten as 



F(x,T) = (<^(x,0)|exp 



6{x,t)dt 







\<P{X,T)) 

(25) 







exp / 5(x,t)dt + 7(x,r) 
Jo 

Hence, when we write the cumulant generating function as niXi't) = log-F(X;^); the 
cumulant generating function is adequately divided into two parts: 

/i(x,T)= / 5(x,t)dt + 7(x,T) 

T 

5(x,t)dt + ln((^(x,0)|</)(x,T)). (26) 



Note that there is no need to calculate the geodesic curve explicitly; the solutions of 
( !T4|) and (fT5|) is enough to calculate all quantities. We will call the first term in (!26|l 
'the dynamical quantity', and the second term 'the geometric quantity'. (Someone feels 
the word 'phase' is strange in this case because the fiber is not U{1) but M_(_, so that we 
here avoid the word 'phase'.) 
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If the time-evolution operator H is time-independent, the ket vector \ip{x, t)) finally 
converges to the eigenstate with the largest eigenvalue, and S{x,t) is characterized by 
the largest eigenvalue; it has been well-known that the cumulant generating function 
is related to the largest eigenvalue of H, and therefore the obtained results from the 
fiber bundle is consistent with the previous works. In addition, t)) = when the 

system is in the eigenstate. This means that the geometric quantity 7(x, T) converges 
to a certain value. Of course, when we consider a nonequilibrium system, there is 
a current which stems from the dynamical quantity, and the dynamical quantity is 
changing continuously with time. The above division shows that there is a certain 'finite' 
contribution from the change of the system, and we can consider the contribution from 
the geometric quantity as an excess quantity discussed a lot in previous works, due to 
the finiteness of the quantity. 



4. Numerical results: Excess quantity 

In order to investigate the characteristics of the geometric quantity, we here perform a 
numerical experiment for a specific model. The following model is considered here, which 
has been used in order to study the geometric contribution under cyclic perturbations 
[9HI2]: 













L 




container 




R 



(27) 



The system consists of three parts. The container can contain at most one particle in 
it, and then the container has only two states, i.e., filled state or empty one. When 
the container is filled with one particle, the particle can escape from the container 
by jumping into one of the two absorbing states: the left reservoir or the right one. In 
contrast, when the container is empty, the left or right reservoirs can emit a new particle 
into the container. The master equation for the system is written as follows: 

— ( P^^^^ \ = ( ~^21-«^21 ^2 + ^?2 ^ ( Plit) \ (28) 
dt \ P2{t) J \ Wl, + Wl, -Wl^ - J \ P2{t) J' 

where pi(t) and P2{t) = l—pi{t) correspond to the probabilities with which the container 
is empty and filled, respectively. 

Here, we calculate the number of hopping from the container to the right reservoir. 
Hence, the quantities {a^s'} is selected as follows: = 1 a'^^, = for other cases. 
The time-evolution operator H{x,t) is therefore written as 

Hix,t)=(-i\i^ -- + --f y (29) 

\^ W2i+W^il -^2-^12 J 

We consider relaxation processes of the system. That is, the transition rates {w'^^,} 
are assumed to be time- independent, and we investigate effects of the initial probability 
Pi = pi(t = 0) . If the probability pi is the stationary one for the master equation fl28|) . 



Excess quantity derived from fiber bundles in counting statistics 



8 



0.003 



= 0.60 
= 0.59 
= 0.63 




0.002 - 



Pi 
Pi 
Pi 




-0.001 



-0.3 -0.2 -0.1 



0.0 
X 



0.1 



0.2 



0.3 



Figure 1. Numerical results for the geometric quantity 7(x,r). Results with 
three different initial probabilities pi(t = 0) are shown, pi ~ 0.60 corresponds to 
the stationary probability of the system. 



the system is staying in the stationary state. On the other hand, if the probabihty pi 
is not stationary one, the system relaxes to the stationary state. The non-stationary 
initial pi corresponds to a stationary state with a different transition matrix, and hence 
the above situation can be considered as the case in which there are sudden changes for 
the transition rates at time t = by external forces. 

Note that the geometric quantity exists even in the stationary initial states. That is, 
the initial state is \p{t = 0)) for (IHj), and (1| for ( fT5l) . and these states are not the right 



and left eigenstates of H{x,t). Hence, |0(x,O)) 7^ \4>{x,T)) and {(p{x,0)\(f){x,T)) 7^ 1 



in general. In the following numerical experiments, we will check whether the geometric 
quantity for the stationary initial cases gives a current or not. 

Setting = 2 and w'^^, = 1 for other cases, we numerically evaluate the geometric 
quantity 7(x, T). Here, the final time T is selected large enough for the relaxation. 
Figure [1] shows the results. When we set the initial probability pi{t = 0) as 0.60, 
which corresponds to the stationary probability of the given transition rates, the first 
derivative at x = is zero. This means that there is no contribution to the current from 
the geometric quantity, and it is reasonable. On the other hand, non-stationary initial 
probabilities pi{t = 0) give a current due to the geometric quantity; the first derivative 
at X = is non-zero. Note that the geometric quantities are calculated for large T, and 
we have confirmed the convergence of these quantities. Hence, it is possible to consider 
the contribution from the geometric quantity as an excess quantity. 

Finally, we comment on the higher-order contribution from the geometric quantity. 
Figure [1] shows that there is a contribution to the current fluctuation of the system 
from the geometric quantity, even when the system starts from the stationary state. 
One may consider that it is strange, because the current and its statistics in the 
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stationary state may be characterized only by the largest eigenvalue of the time-evolution 
operator H{x,t). However, it has been known that there is another source of the 
contribution for the current and its statistics; for example, so-called 'boundary terms' 
appear when we formulate the counting statistics using the path-integrals. There are at 
least two different schemes for the boundary terms p^lTG] : in one scheme, the boundary 
terms are included in the dynamical quantity [15]; in another scheme, they are in the 
geometric quantity [16] . The formulation based on the fiber bundles suggests a different 
consequence; the boundary terms are divided into two parts in some way and contribute 
to both the dynamical and geometric quantities. 

5. Concluding remarks 

In the present paper, we proposed a natural division of the cumulant generating function 
for the current. The discussion is based on the fiber bundles and, using a numerical 
experiment, we confirmed that the geometric quantity introduced by the fiber bundle 
discussions has a characteristics of excess quantities. 

We comment on some unclear points. In the numerical experiment, it has been 
revealed that even in the stationary state the geometric quantity has the contribution to 
the fluctuations. The physical meaning of this fact should be investigated in future works 
although the contribution from the geometric quantity will be very small comparing to 
that from the dynamical quantity in general. In addition, a derivation of the previous 
path-integral results for adiabatic cases [T5|[T6] from the fiber bundle formulation is still 
unclear. Mathematically more rigorous and extended discussions may be also needed; 
the non-Hermite property could cause difficulty for the settings of the fiber bundle and 
its dual. A cyclic time-evolution for non-Hermitian cases has been discussed in [23] . 
and the settings in the present paper seem to adequately correspond to the formulation 
in [23]. As far as we checked, the mathematical discussion in the present paper is valid 
and adequate when we consider only the actual time-evolution and the geodesic curves. 
However, there may be more adequate discussions for the non-Hermite cases, which 
enables us to give extended discussions. 
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Appendix A. Comment on the geodesic curve 

Consider the trajectory in P, which is made by the time-evolution from the initial ket 
vector t = 0)) to the final ket vector t = T)). In order to make a closed curve 
in P, a concept of 'geodesic' curves is useful [22] . 
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Firstly, we explain how to derive the geodesic curve. Considering a curve in P, 
s)), we introduce the tangent vector to this curve as |n) = ■j^\(t>{x^ s))^ discussed 
in section 2. In addition, we also introduce the corresponding tangent vector in the dual 
space as {v\ = -^{ip{x, s)\. In the following discussions, we will use A'^{x,s) in order to 
express the connection in the dual space, but we assume that the dual curve is selected 
so that the inner produce is conserved; i.e., A{x, s) = s) as discussed in section 2. 

Under a gauge transformation, the tangent vectors \u) and {v\ are not gauge covariant, 
and therefore the following gauge covariant quantities are useful: 

\u'{x, s)) ^ s)) - A{x, sMix. s)), (A.l) 

{v'ix,s)\ ^ ^Jv^ix,s)\ + {v^ix.s)\A\x,s), (A.2) 

Hence, {v'\u') is gauge invariant, and {v'\u') is available in order to define a kind of 
metric. The variation of J {v'\u')dl, where / is an affine parameter, gives the geodesic 
curve s)) and {(p{x^s)\, which are determined by 

±\u'{x,s))-A^{x,s)\u'{x,s)) = 0, (A.3) 

±{v'{x,s)\ + {v'{x,s)\A{x,s)=0. (A.4) 

Using the geodesic curve, it is possible to evaluate the integral ( IMl) . We consider a 
geodesic curve \4>^{x, q)) {q ^ [0, 1]) connecting to the trajectory obtained from the time- 
evolution backward; \(p^ix,q = 0)) = |0(x,^ = T)) and |0^(x,g = 1)) = |0(x,^ = 0)). 
Hence, we have a closed curve in P, i.e., the actual time-evolution from \4>{x,t = 0)) 
to \(p{Xyt = T)), and back along the geodesic curve \(p^{x^Q))- We can therefore define 
the integral ^ A{x,s)ds along the closed curve. Note that A{x,t) = on the actual 
time-evolution, and A{x, s) 7^ on \(f)^{x, q)) in general. In order to calculate the value 
of the integral, the characteristics of the geodesic curve is important. Using a certain 
gauge transformation on \(t>^{x^Q)) = 6xp(— /3(g))|0§(x, g)), where /3{q = 0) = 0, we 
choose a specific connection A{x, g) = on |0^(x, q))- Due to the gauge transformation, 
the derived whole curve, including the actual time-evolution, is not closed. Note that 
flA.Sp and flA.4l) are gauge covariant, and the gauge transformed curve \(f>^{XyQ)) is 
still a geodesic curve. In order to have a closed curve, we must use a vertical curve 
joining \4>^{x^Q = 1)) and \4>{Xit = 0)). It is clear that the vertical contribution 
corresponds to the integral § A{x, s)ds, because A{x, t) = on the actual time-evolution 
and A{x, s) = on the gauge transformed curve. In addition, the contribution from the 
vertical curve is given by f3{q = 1). In order to evaluate f3{q = 1), the quantity g{q) = 
((^^(x, 5')|0(X! ^ = T)), where (<^^(X)^)I is the dual of \(p^{XyQ))y is considered. Then, 
we have g{0) = 1 because {(p^{xA = O)|0(x,^ = T)) = {(p^{x,q = O)|0(x,i = T)) = 1; 
■^g{0) = 0, which is obtained from A{x, ? = 0) = 0; and -^g{q) = 0, which stems from 
the definitions of the geodesic curve. These facts suggest that g{q) is constant for all q. 
Hence, 

1 = {^'ix, q = 0) |0(X, t = T)) = mx, q = 0) |0(x, t = T)) 
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= e-^^'^='H^ix,t = OMix,t = T)). (A.5) 

That is, the vertical contribution is evaluated as (3{q = 1) = \n{(p{x, t = 0)\(f){xy t = T)), 
which gives ( l24l) . 

References 

Oono Y and Paniconi M 1998 Prog. Theor. Phys. Supple. 130 29 
Hatano T and Sasa S 2001 Phys. Rev. Lett. 86 3463 
Speck T and Seifert U 2005 J. Phys. A: Math. Gen. 38 L581 
Esposito M, Harbola U and Mukamel S 2007 Phys. Rev. E 76 031132 
Komatsu T S and Nakagawa N 2008 Phys. Rev. Lett. 100 030601 
Komatsu T S, Nakagawa N, Sasa S I and Tasaki H 2008 Phys. Rev. Lett. 100 230602 
Komatsu T S, Nakagawa N, Sasa S I and Tasaki H 2009 J. Stat. Phys. 134 401 
Esposito M and den Broeck C V 2010 Phys. Rev. Lett. 104 090601 
Sinitsyn N A and Nemenman I 2007 Europhys. Lett. 77 58001 
Sinitsyn N A and Nemenman I 2007 Phys. Rev. Lett. 99 220408 
Ohkubo J 2008 J. Stat. Mech. P02011 
Ohkubo J 2008 J. Chem. Phys. 129 205102 

Sinitsyn N A and Saxena A 2008 J. Phys. A: Math. Theor. 41 392002 
Sinitsyn N A 2009 J. Phys. A: Math. Theor. 42 193001 
Sinitisyn N A and Nemenman I 2010 lET Syst. Biol. 4 409 
Sagawa T and Hayakawa H 2011 Phys. Rev. E 84 051110 
Ohkubo J and Eggel T 2010 J. Phys. A: Math. Theor. 43 425001 
Gopich I V and Szabo A 2006 J. Chem. Phys. 124 154712 

Garrahan J P, Jack R L, Lecomte V, Pitard E, van Duijvcndijk K and van Wijland F 2009 J. 

Phys. A: Math. Theor. 42 075007 
Bohm A, Mostafazadeh A, Koizumi H, Niu Q and Zwanziger J 2003 The Geometric Phase in 

Quantum Systems (Berlin: Springer- Verlag) 
Nakahara M 2003 Geometry, Topology, and Physics 2nd ed. (Boca Raton: Taylor & Francis) 
Samuel J and Bhandari R 1988 Phys. Rev. Lett. 60 2339 
Garrison J C and Wright E M 1988 Phys. Lett. A 128 177 



